We want to characterize the CNV data from St. Jude OS samples in order to determine if the CNVs are consistent across samples within a patient. This will tell us if these tumors are genomically stable or if there are ongoing rearrangements. We also want to determine if there are sub-clones within a single patient which would be revealed by distinct CNV patterns within a patient.
For this analysis we got access to many St. Jude OS samples. We also included our single-cell sequencing data which we treated as bulk for the purposes of these analyses.
library(tidyverse)
## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'
## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'hms'
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.3.4 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.7
## ✔ tidyr 1.1.3 ✔ stringr 1.4.0
## ✔ readr 2.0.1 ✔ forcats 0.5.1
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
theme_set(theme_bw())
theme_update(plot.title = element_text(hjust = 0.5))
library(gt)
for directoryName in \
misc \
slurmOut \
input\split \
output \
output/aligned \
output/figures \
output/counts \
output/varscan/copynumber/ \
output/vcfs \
output/vcfs/realignSummary \
output/vcfs/p53Char \
output/vcfs/p53Char/clinVar
do
if [ ! -d ${directoryName} ]
then
mkdir -p ${directoryName}
fi
done
Had to do this by hand due to the links expiring after 24 hours
We had a lot of trouble downloading large bam files from DNAnexus. The downloads frequently got interupted. To get around this, I split the files in DNAnexus into smaller chunks to make interuptions less likely. I also wrote the code below to automatically resume the downloads if the download was interupted. I had the script check the md5sums to be sure that the file wasn’t corrupted during download. The original file is then re-generated using cat and md5sum checked again.
Don’t forget to dos2unix the .txt files before running!
perl ../scripts/prepDownloadFile.pl \
--linkFile misc/DNAnexus_export_urls-20210926-103118.txt \
--md5File misc/allMd5.txt \
> misc/batch1DownloadFile.txt
perl ../scripts/prepDownloadFile.pl \
--linkFile misc/DNAnexus_export_urls-20210927-085055.txt \
--md5File misc/allMd5.txt \
> misc/batch2DownloadFile.txt
perl ../scripts/prepDownloadFile.pl \
--linkFile misc/DNAnexus_export_urls-20210928-082401.txt \
--md5File misc/allMd5.txt \
> misc/batch3DownloadFile.txt
perl ../scripts/prepDownloadFile.pl \
--linkFile misc/DNAnexus_export_urls-20210929-083226.txt \
--md5File misc/allMd5.txt \
> misc/batch4DownloadFile.txt
perl ../scripts/prepDownloadFile.pl \
--linkFile misc/DNAnexus_export_urls-20210930-084226.txt \
--md5File misc/allMd5.txt \
> misc/batch5DownloadFile.txt
perl ../scripts/prepDownloadFile.pl \
--linkFile misc/DNAnexus_export_urls-20211001-084849.txt \
--md5File misc/allMd5.txt \
> misc/batch6DownloadFile.txt
I had to run these one by one manually due to how the data are organized
sbatch sbatchCmds/forceDownload_1.sh
sbatch sbatchCmds/forceDownload_2.sh
sbatch sbatchCmds/forceDownload_3.sh
sbatch sbatchCmds/forceDownload_4.sh
sbatch sbatchCmds/forceDownload_5.sh
sbatch sbatchCmds/forceDownload_6.sh
cat misc/rawMd5s/rawBamMd5_md5sum_job-G5FQ* > misc/rawBamMd5s.txt
sbatch sbatchCmds/catSplitFiles.sh
for original in S0113 S0114 S0115 S0116 S0126 S0127 S0128
do
ln -s \
/home/gdrobertslab/lab/Counts/${original}/possorted_bam.bam \
input/${original}.bam
done
# Some of the samples were already in the input folder and don't need symLinks
for original in \
SJOS046149_R1 SJOS046149_R2 SJOS046149_X1 SJOS031478_D1 \
SJOS031478_D2 SJOS031478_D3
do
ln -s \
~/analyses/roberts/stjude/rawData/${original}.WholeGenome.bam \
input/${original}.WholeGenome.bam
done
Kicked out the unknown and random chromosomes
cat /reference/homo_sapiens/hg38/ucsc_assembly/illumina_download/Sequence/WholeGenomeFasta/genome.fa.fai | \
grep -v "chrUn\|random\|EBV" \
> misc/chrList.txt
if [ ! -d /gpfs0/scratch/mvc002/roberts/stjude/realigned ]
then
mkdir -p /gpfs0/scratch/mvc002/roberts/stjude/realigned
fi
sbatch sbatchCmds/realign.sh
sbatch sbatchCmds/mpileupVarscan.sh
sbatch sbatchCmds/varScanCopycaller.sh
pdf(file = "output/figures/CNVdistribution_copynumber.pdf")
for (sample_name in list.dirs(path = "output/varscan/copynumber",
full.names = FALSE,
recursive = FALSE)) {
print(read_tsv(list.files(path = paste("output/varscan/copynumber/",
sample_name,
sep = ""),
pattern = "*.copynumber",
full.names = TRUE),
show_col_types = FALSE) %>%
sample_n(1000000) %>%
select(log2_ratio) %>%
mutate(ratio = 2^log2_ratio) %>%
ggplot(aes(x = ratio)) +
geom_histogram(bins = 300) +
geom_vline(xintercept = 1) +
scale_x_log10() +
geom_vline(xintercept = 2) +
geom_vline(xintercept = 0.5) +
ggtitle(sample_name))
}
dev.off()
## png
## 2
This approach has the issue that normalization isn’t perfect across the samples due to the extensive number of genomic changes. This artificially increases the distance between similar samples. This likely will not be used in the paper.
perl scripts/varscanToMatrix.pl \
-v \
--fileList "output/varscan/copycaller/S*/*chr*.txt" \
> output/varscan/varscan_distMatrix.txt
perl scripts/varscanToMatrix_threshhold.pl \
-v \
--fileList "output/varscan/copycaller/S011*/*chr*.txt" \
> output/varscan/varscan_distMatrix_threshhold.txt
sbatch sbatchCmds/combineVarscan.sh
large_bin_size <- 100000
combined_data <- read_tsv("output/varscan/combined_calls.txt",
show_col_types = FALSE) %>%
mutate(chr_num = str_remove(chr, "chr") %>%
as.numeric()) %>%
arrange(chr_num, bin) %>%
mutate(order = seq_len(nrow(.)),
larger_bin = floor(bin / large_bin_size) * large_bin_size) %>%
pivot_longer(cols = c(-chr, -bin, -order, -chr_num, -larger_bin),
names_to = "sample",
values_to = "log_ratio")
patient_key <- data.frame(
patient = unique(combined_data$sample) %>%
str_remove("_.+") %>%
str_replace("S0126", "SJOS030645") %>%
str_replace("S0127", "SJOS030645") %>%
str_replace("S0128", "SJOS030645") %>%
str_replace("S0114", "SJOS031478") %>%
str_replace("S0116", "SJOS031478") %>%
str_replace("S0113", "SJOS046149") %>%
str_replace("S0115", "SJOS046149"),
sample_type = unique(combined_data$sample) %>%
str_replace("^S0.+", "Xenograft") %>%
str_replace(".+_X.+", "Xenograft") %>%
str_replace(".+_D.+", "Diagnosis") %>%
str_replace(".+_M.+", "Metastasis") %>%
str_replace(".+_R.+", "Relapse"),
sample = unique(combined_data$sample))
combined_data <- full_join(combined_data, patient_key)
## Joining, by = "sample"
chr_cols <- sample(rainbow(length(unique(combined_data$chr))),
length(unique(combined_data$chr)),
replace = FALSE)
for (patient_name in unique(patient_key$patient)) {
plot_name <- combined_data %>%
filter(patient == patient_name) %>%
group_by(chr, sample, larger_bin) %>%
summarize(log_ratio = median(log_ratio),
order = min(order),
.groups = "drop") %>%
ggplot(aes(x = order,
y = log_ratio,
color = chr)) +
geom_hline(yintercept = 0, color = "black") +
geom_point(alpha = 0.5,
size = 0.5) +
scale_color_manual(values = chr_cols) +
labs(y = "Copy number log ratio",
x = "") +
facet_wrap(~ sample, ncol = 1)
ggsave(paste0("output/figures/varscan",
patient_name,
".png"),
plot = plot_name,
width = 20,
height = 8)
}
## Take a look at regions commonly amplified/deleted across samples
unique_clones <- c("SJOS001101_M4",
"SJOS001111_M1",
"SJOS001116_X2",
"SJOS001121_X2",
"SJOS001126_X1",
"SJOS013768_X1",
"SJOS016016_X1",
"SJOS030101_R1",
"SJOS030645_D2",
"SJOS031478_D3",
"SJOS046149_R2",
"SJOS001105_R1")
combined_data %>%
filter(sample %in% unique_clones) %>%
group_by(chr, bin) %>%
summarize(mean_ratio = median(log_ratio),
order = min(order),
.groups = "drop") %>%
ggplot(aes(x = bin,
y = mean_ratio,
color = chr)) +
geom_hline(yintercept = 0, color = "black") +
geom_point(alpha = 0.5,
size = 0.5) +
labs(x = "",
y = "Mean copy number ratio") +
facet_wrap(~ chr, scales = "free_x") +
theme(legend.position = "none")
ggsave("output/figures/clonal_mean_CNV.png",
width = 15,
height = 10)
high_chr3 <- combined_data %>%
filter(sample %in% unique_clones) %>%
group_by(chr, bin) %>%
summarize(mean_ratio = median(log_ratio),
order = min(order),
.groups = "drop") %>%
filter(chr == "chr3" & mean_ratio < -0.5)
cnv_cor <- read_tsv("output/varscan/combined_calls.txt",
col_names = TRUE,
show_col_types = FALSE) %>%
mutate(chr_bin = paste(chr, bin, sep = "_"), .keep = "unused") %>%
column_to_rownames("chr_bin") %>%
cor()
patient_key <- data.frame(patient = rownames(cnv_cor) %>%
str_remove("_.+") %>%
str_replace("S0126", "SJOS030645") %>%
str_replace("S0127", "SJOS030645") %>%
str_replace("S0128", "SJOS030645") %>%
str_replace("S0114", "SJOS031478") %>%
str_replace("S0116", "SJOS031478") %>%
str_replace("S0113", "SJOS046149") %>%
str_replace("S0115", "SJOS046149"),
sample_type = rownames(cnv_cor) %>%
str_replace("^S0.+", "Xenograft") %>%
str_replace(".+_X.+", "Xenograft") %>%
str_replace(".+_D.+", "Diagnosis") %>%
str_replace(".+_M.+", "Metastasis") %>%
str_replace(".+_R.+", "Relapse"),
sample_2 = rownames(cnv_cor)) %>%
column_to_rownames("sample_2")
long_cor <-
cnv_cor %>%
as.data.frame() %>%
rownames_to_column("sample_1") %>%
pivot_longer(-sample_1,
names_to = "sample_2",
values_to = "cor") %>%
left_join(patient_key %>%
rownames_to_column("sample_1") %>%
rename(patient_1 = patient) %>%
select(-sample_type)) %>%
left_join(patient_key %>%
rownames_to_column("sample_2") %>%
rename(patient_2 = patient) %>%
select(-sample_type)) %>%
mutate(within = patient_1 == patient_2,
sample_1 = fct_reorder(sample_1, patient_1),
sample_2 = fct_reorder(sample_2, patient_2)) %>%
filter(sample_1 != sample_2) %>%
arrange(patient_1)
## Joining, by = "sample_1"Joining, by = "sample_2"
## Try to id clones
clone_table <- tibble(sample_1 = character(),
clone_num = factor())
claimed_list <- character()
cor_cutoff <- 0.4
clone_num <- 1
for (sample_name in rownames(cnv_cor)) {
if (!sample_name %in% claimed_list) {
new_claimed <-
long_cor %>%
filter(sample_1 == sample_name &
cor >= cor_cutoff &
within == TRUE) %>%
select(sample_2, cor) %>%
mutate(clone_num = as.factor(clone_num)) %>%
bind_rows(tibble(clone_num = as.factor(clone_num),
sample_2 = sample_name,
cor = 1))
clone_table <-
rbind(clone_table, new_claimed)
claimed_list <- c(claimed_list, sample_name, new_claimed$sample_2)
clone_num <- clone_num + 1
}
}
clone_table <-
clone_table %>%
group_by(sample_2) %>%
arrange(cor * -1, .by_group = TRUE) %>%
slice_head(n = 1) %>%
ungroup() %>%
select(-cor)
long_cor <- full_join(long_cor, clone_table, by = "sample_2")
# Mean within patient correlation
long_cor %>%
filter(within == TRUE) %>%
pull(cor) %>%
summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.007455 0.362223 0.509682 0.503482 0.717656 0.921232
# # Mean correlation within a clone
# long_cor %>%
# select(sample_2, patient_2, cor, clone_num) %>%
# group_by(patient_2, clone_num) %>%
# summarize(cor = mean(cor)) %>%
# pull(cor) %>%
# summary()
# Get number of patients with >1 clone
long_cor %>%
select(patient_2, clone_num) %>%
arrange(patient_2) %>%
distinct() %>%
pull(patient_2) %>%
table()
## .
## SJOS001101 SJOS001105 SJOS001111 SJOS001116 SJOS001121 SJOS001126 SJOS013768
## 1 1 1 1 1 1 2
## SJOS016016 SJOS030101 SJOS030645 SJOS031478 SJOS046149
## 1 2 2 2 2
# Mean within clone correlation
long_cor %>%
full_join(clone_table %>%
rename(sample_1 = sample_2,
clone_2 = clone_num)) %>%
filter(clone_num == clone_2) %>%
pull(cor) %>%
summary()
## Joining, by = "sample_1"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4007 0.5312 0.6807 0.6707 0.7850 0.9212
# Mean correlatin within a patient between divergent clones
long_cor %>%
full_join(clone_table %>%
rename(sample_1 = sample_2,
clone_2 = clone_num)) %>%
filter(patient_1 == patient_2 &
clone_num != clone_2) %>%
pull(cor) %>%
summary()
## Joining, by = "sample_1"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.007455 0.123653 0.307508 0.282802 0.409872 0.584116
# Mean correlation between patients
long_cor %>%
full_join(clone_table %>%
rename(sample_1 = sample_2,
clone_2 = clone_num)) %>%
filter(patient_1 != patient_2) %>%
pull(cor) %>%
summary()
## Joining, by = "sample_1"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.1363 0.1165 0.1823 0.1781 0.2438 0.4527
# Correlation between xenografts and original samples within a single clone and patient
long_cor %>%
full_join(clone_table %>%
rename(sample_1 = sample_2,
clone_2 = clone_num)) %>%
left_join(patient_key %>%
rownames_to_column("sample_1") %>%
rename(type_1 = sample_type) %>%
select(-patient)) %>%
left_join(patient_key %>%
rownames_to_column("sample_2") %>%
rename(type_2 = sample_type) %>%
select(-patient)) %>%
filter(patient_1 == patient_2 &
clone_num == clone_2 &
type_1 == "Xenograft" &
type_2 != "Xenograft") %>%
select(-within, -patient_1, -patient_2) %>%
pull(cor) %>%
summary()
## Joining, by = "sample_1"Joining, by = "sample_1"Joining, by = "sample_2"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4710 0.5622 0.7224 0.7019 0.8295 0.9212
ggplot(long_cor, aes(x = cor)) +
geom_histogram(bins = 200) +
facet_wrap(~ within, ncol = 1)
# correlation of about 0.4 seems to be the upper bound of non-related samples
n_clones <- max(as.numeric(long_cor$clone_num))
text_col <-
long_cor %>%
full_join(tibble(clone_num = seq_len(n_clones) %>%
as.factor(),
text_col = rainbow(n = n_clones) %>%
sample(size = n_clones))) %>%
mutate(group = paste(patient_2, sample_2)) %>%
select(group, text_col) %>%
distinct() %>%
arrange(group) %>%
pull(text_col)
## Joining, by = "clone_num"
ggplot(long_cor, aes(x = paste(patient_1, sample_1),
y = paste(patient_2, sample_2),
fill = cor)) +
geom_tile() +
scale_fill_gradient2(low = "blue",
mid = "white",
high = "red",
midpoint = 0.5) +
theme(axis.text.y = element_text(colour = text_col))
## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
patient_key <- patient_key %>%
rownames_to_column("sample_2") %>%
full_join(clone_table) %>%
column_to_rownames("sample_2") %>%
rename(clone = clone_num,
`sample type` = sample_type) %>%
relocate(`sample type`, .after = last_col())# %>%
## Joining, by = "sample_2"
# select(-`sample type`)
col_patient <- list(patient = c("#a6cee3",
"#1f78b4",
"#b2df8a",
"#33a02c",
"#fb9a99",
"#e31a1c",
"#fdbf6f",
"#ff7f00",
"#cab2d6",
"#6a3d9a",
"#ffff99",
"#b15928"),
clone = c("#96FF00",
"#00FFD2",
"#FF00B4",
"#2b754e",
"#3CFF00",
"#FF5A00",
"#F000FF",
"#9600FF",
"#3C00FF",
"#F0FF00",
"#FFB400",
"#0078FF",
"#001EFF",
"#00FF1E",
"#FF0000",
"#00D2FF",
"#FF005A"),
`sample type` = c("#a6cee3",
"#1f78b4",
"#b2df8a",
"#33a02c"))
names(col_patient$patient) <- patient_key$patient %>%
as.factor() %>%
levels()
names(col_patient$`sample type`) <- patient_key$`sample type` %>%
as_factor() %>%
levels()
names(col_patient$clone) <- levels(clone_table$clone_num)
png("output/figures/cnv_correlation_heatmap.png",
width = 3500,
height = 3000,
res = 300)
pheatmap::pheatmap(cnv_cor,
annotation_row = patient_key,
annotation_col = patient_key,
annotation_colors = col_patient)
dev.off()
## png
## 3
Filter out indels and SNPs with qual < 20 and depth < 20
# Make up a list of chromosomes to analyze for use by parallel
grep "^>" \
/reference/homo_sapiens/hg38/ucsc_assembly/illumina_download/Sequence/WholeGenomeFasta/genome.fa \
| grep -v "random\|chrUn\|chrEBV" \
| perl -pe 's/>//' \
> misc/chrList.txt
sbatch sbatchCmds/callVariants.sh
sbatch sbatchCmds/indexVcfs.sh
module load GCC/9.3.0 \
GCCcore/9.3.0 \
BCFtools/1.11 \
SAMtools/1.15
perl scripts/vcfToMatrix_parallel.pl \
--fileList "output/vcfs/S*vcf.gz" \
--threads 20 \
> output/vcfs/snpDistanceMatrix.txt
sbatch sbatchCmds/TP53MutChar.sh
wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar_20220416.vcf.gz
mv clinvar_20220416.vcf.gz misc/
zcat misc/clinvar_20220416.vcf.gz \
| awk '$1 ~ /^#/ || \
$1 == 17 && \
$2 >= 7600000 && \
$2 <= 7690000 \
{print}' \
| perl -pe 's/^17/chr17/' \
| gzip \
> misc/clinvar_TP53.vcf.gz
combineClinvar.sh
filter_list <- function(x) {
str_split(x, ",") %>%
unlist() %>%
str_subset("TP53") %>%
str_c(collapse = ",")
}
file_list <- list.files("output/vcfs/p53Char/clinVar",
pattern = ".+_G.+",
full.names = TRUE)
combined_data <- tibble(sample = character())
for (file_name in file_list) {
combined_data <-
read_tsv(file_name,
col_names = TRUE,
col_types = cols(.default = "c")) %>%
filter(grepl("TP53", ANN)) %>%
rowwise() %>%
mutate(ANN = filter_list(ANN),
sample = str_match(file_name, ".+/(.+).txt")[2]) %>%
filter(ANN != "") %>%
bind_rows(combined_data) %>%
relocate(sample, CLNSIG, ANN)
}
write_tsv(combined_data,
file = "output/vcfs/p53Char/p53Char_combined_germline.txt")
perl scripts/cmpVarscanCalls.pl \
--binSize 100 \
--inputFiles "output/varscan/copycaller/*/*chr17.txt" \
> output/varscan/combined_calls_chr17.txt
##
## Printing out results
P53_loc <- list(min = 7661779,
max = 7687538)
combined_data <- read_tsv("output/varscan/combined_calls_chr17.txt",
show_col_types = FALSE) %>%
mutate(chr_num = str_remove(chr, "chr") %>%
as.numeric()) %>%
filter(chr_num == 17 & bin >= 7550000 & bin <= 7800000) %>%
arrange(chr_num, bin) %>%
pivot_longer(cols = c(-chr, -bin, -chr_num),
names_to = "sample",
values_to = "log_ratio")
patient_key <- data.frame(
patient = unique(combined_data$sample) %>%
str_remove("_.+") %>%
str_replace("S0126", "SJOS030645") %>%
str_replace("S0127", "SJOS030645") %>%
str_replace("S0128", "SJOS030645") %>%
str_replace("S0114", "SJOS031478") %>%
str_replace("S0116", "SJOS031478") %>%
str_replace("S0113", "SJOS046149") %>%
str_replace("S0115", "SJOS046149"),
sample_type = unique(combined_data$sample) %>%
str_replace("^S0.+", "Xenograft") %>%
str_replace(".+_X.+", "Xenograft") %>%
str_replace(".+_D.+", "Diagnosis") %>%
str_replace(".+_M.+", "Metastasis") %>%
str_replace(".+_R.+", "Relapse"),
sample = unique(combined_data$sample))
combined_data <- full_join(combined_data, patient_key) %>%
mutate(amplification = if_else(log_ratio > 0.5,
"amplification",
if_else(log_ratio < -0.5,
"deletion",
"normal")))
## Joining, by = "sample"
ggplot(combined_data,
aes(x = bin,
y = sample,
fill = log_ratio)) +
geom_tile() +
scale_fill_gradient2(low = "#8e0152",
mid = "white",
high = "#276419") +
facet_grid(patient ~ .,
scales = "free_y",
space = "free") +
theme(strip.text.y.right = element_text(angle = 0)) +
geom_vline(xintercept = P53_loc$min,
linetype = "dashed",
colour = "red") +
geom_vline(xintercept = P53_loc$max,
linetype = "dashed",
colour = "red") +
labs(x = "Chromosome 17 position",
y = "",
title = "TP53 copy numbers")
ggsave(file = "output/figures/p53Cnv_chr17.png",
width = 10,
height = 10)
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux 8.1 (Ootpa)
##
## Matrix products: default
## BLAS: /gpfs0/export/apps/opt/R/4.1.0-foss-2020a/lib64/R/lib/libRblas.so
## LAPACK: /gpfs0/export/apps/opt/R/4.1.0-foss-2020a/lib64/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gt_0.3.0 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7
## [5] purrr_0.3.4 readr_2.0.1 tidyr_1.1.3 tibble_3.1.7
## [9] ggplot2_3.3.4 tidyverse_1.3.1 rmarkdown_2.9
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.1 xfun_0.24 bslib_0.2.5.1 haven_2.4.1
## [5] colorspace_2.0-1 vctrs_0.3.8 generics_0.1.0 htmltools_0.5.1.1
## [9] yaml_2.2.1 utf8_1.2.1 rlang_1.0.2 jquerylib_0.1.4
## [13] pillar_1.7.0 glue_1.4.2 withr_2.4.2 DBI_1.1.1
## [17] dbplyr_2.1.1 modelr_0.1.8 readxl_1.3.1 lifecycle_1.0.0
## [21] munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0 rvest_1.0.0
## [25] evaluate_0.14 knitr_1.33 tzdb_0.1.2 ps_1.6.0
## [29] fansi_0.5.0 broom_0.7.7 Rcpp_1.0.7 backports_1.2.1
## [33] scales_1.1.1 jsonlite_1.7.2 fs_1.5.0 hms_1.1.0
## [37] digest_0.6.27 stringi_1.6.2 grid_4.1.0 cli_2.5.0
## [41] tools_4.1.0 magrittr_2.0.1 sass_0.4.0 crayon_1.4.1
## [45] pkgconfig_2.0.3 ellipsis_0.3.2 xml2_1.3.2 reprex_2.0.0
## [49] lubridate_1.7.10 rstudioapi_0.13 assertthat_0.2.1 httr_1.4.2
## [53] R6_2.5.0 compiler_4.1.0